source('00_util_scripts/mod_bulk.R')

late.kc <- c('Krt1','Krt10','Lor','Ivl','Tgm1','Flg') |>
  str_to_upper()

# mva pathway ----------
# ACAT >> HMGCS >> HMGCR >> MVK >> PMVK >> MVD >> IDI >> (FDPS,GGPS1)
kegg_mva <-
  c('Acat1','Acat2','Hmgcs1','Hmgcs2','Hmgcr','Mvk','Pmvk','Mvd','Idi1','Idi2','Fdps') |>
  str_to_upper()

mva_fct <- tibble(gene = kegg_mva, ordered = fct_inorder(kegg_mva))

# interested cytokines ------
key_cytokine <- c('Ccl20','Tslp','Flt3lg','Csf2','Tnf','Il6') |>
  str_to_upper()

# choi24 ---------
choi24 <- download_ncbi_counts('GSE250127')

tdb.choi24 <- choi24 |>
  select(Symbol, GSM7974261, GSM7974262, GSM7974263) |>
  pivot_longer(-1) |>
  mutate(group = ifelse(str_ends(name, '261'), 'ctrl', 'ros')) |>
  tidybulk(.sample = name, .transcript = Symbol, .abundance = value)

tdb.choi24 <- tdb.choi24 |>
  preproc_bulk(group = group)

tdb.choi24 |>
  plot_qc_bulk(scaled_abundance = value_scaled, group)

res.choi24 <- tdb.choi24 |>
  test_differential_abundance(~ 0 + group, omit_contrast_in_colnames = T,
                              contrasts = 'groupros-groupctrl')

res.choi24 |>
  pivot_transcript() |>
  filter(Symbol %in% kegg_mva) |>
  filter(FDR < .1)

tdb.choi24 |>
  filter(Symbol %in% kegg_mva, !str_detect(name, '63$')) |>
  ggplot(aes(group, value_scaled)) +
  geom_col() +
  facet_wrap(~Symbol, scales = 'free_y')
